

<!DOCTYPE html>
<!--[if IE 8]><html class="no-js lt-ie9" lang="en" > <![endif]-->
<!--[if gt IE 8]><!--> <html class="no-js" lang="en" > <!--<![endif]-->
<head>
  <meta charset="utf-8">
  
  <meta name="viewport" content="width=device-width, initial-scale=1.0">
  
  <title>Solutions &mdash; quippy 7ea04c5+ documentation</title>
  

  
  
    <link rel="shortcut icon" href="../_static/favicon.ico"/>
  
  
  

  
  <script type="text/javascript" src="../_static/js/modernizr.min.js"></script>
  
    
      <script type="text/javascript" id="documentation_options" data-url_root="../" src="../_static/documentation_options.js"></script>
        <script type="text/javascript" src="../_static/jquery.js"></script>
        <script type="text/javascript" src="../_static/underscore.js"></script>
        <script type="text/javascript" src="../_static/doctools.js"></script>
        <script type="text/javascript" src="../_static/language_data.js"></script>
        <script async="async" type="text/javascript" src="https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.5/latest.js?config=TeX-AMS-MML_HTMLorMML"></script>
        <script type="text/x-mathjax-config">MathJax.Hub.Config({"tex2jax": {"processClass": "math|output_area", "processEscapes": true, "ignoreClass": "document", "inlineMath": [["$", "$"], ["\\(", "\\)"]]}})</script>
    
    <script type="text/javascript" src="../_static/js/theme.js"></script>

    

  
  <link rel="stylesheet" href="../_static/css/theme.css" type="text/css" />
  <link rel="stylesheet" href="../_static/pygments.css" type="text/css" />
  <link rel="stylesheet" href="../_static/theme_overrides.css" type="text/css" />
    <link rel="index" title="Index" href="../genindex.html" />
    <link rel="search" title="Search" href="../search.html" />
    <link rel="next" title="References" href="adaptive-qmmm-references.html" />
    <link rel="prev" title="Step 3: LOTF hybrid MD simulation of fracture in Si" href="adaptive-qmmm-step3.html" /> 
</head>

<body class="wy-body-for-nav">

   
  <div class="wy-grid-for-nav">
    
    <nav data-toggle="wy-nav-shift" class="wy-nav-side">
      <div class="wy-side-scroll">
        <div class="wy-side-nav-search" >
          

          
            <a href="../index.html" class="icon icon-home"> quippy
          

          
            
            <img src="../_static/hybrid.png" class="logo" alt="Logo"/>
          
          </a>

          
            
            
              <div class="version">
                7ea04c5+
              </div>
            
          

          
<div role="search">
  <form id="rtd-search-form" class="wy-form" action="../search.html" method="get">
    <input type="text" name="q" placeholder="Search docs" />
    <input type="hidden" name="check_keywords" value="yes" />
    <input type="hidden" name="area" value="default" />
  </form>
</div>

          
        </div>

        <div class="wy-menu wy-menu-vertical" data-spy="affix" role="navigation" aria-label="main navigation">
          
    
            
            
              
            
            
              <ul class="current">
<li class="toctree-l1"><a class="reference internal" href="../intro.html">Introduction to <code class="docutils literal notranslate"><span class="pre">quippy</span></code></a></li>
<li class="toctree-l1"><a class="reference internal" href="../install.html">Installation of QUIP and quippy</a></li>
<li class="toctree-l1 current"><a class="reference internal" href="index.html">Tutorials</a><ul class="current">
<li class="toctree-l2"><a class="reference internal" href="Introduction.html">Getting started with <code class="docutils literal notranslate"><span class="pre">Atoms</span></code></a></li>
<li class="toctree-l2"><a class="reference internal" href="quippy-ase-interoperability.html">Interoperability with Atomic Simulation Environment</a></li>
<li class="toctree-l2"><a class="reference internal" href="tutorial.html">Introductory Tutorial</a></li>
<li class="toctree-l2"><a class="reference internal" href="quippy-descriptor-tutorial.html">Using Quippy to compute descriptors of atomic environments</a></li>
<li class="toctree-l2"><a class="reference internal" href="adglass.html">Molecular Dynamics Simulation of Fracture in Quartz</a></li>
<li class="toctree-l2 current"><a class="reference internal" href="adaptive-qmmm.html">Adaptive QM/MM MD of Fracture in Silicon</a><ul class="current">
<li class="toctree-l3"><a class="reference internal" href="adaptive-qmmm-step0.html">Introduction</a></li>
<li class="toctree-l3"><a class="reference internal" href="adaptive-qmmm-theory.html">Theoretical background</a></li>
<li class="toctree-l3"><a class="reference internal" href="adaptive-qmmm-step1.html">Step 1: Setup of the Silicon model system</a></li>
<li class="toctree-l3"><a class="reference internal" href="adaptive-qmmm-step2.html">Step 2: Classical MD simulation of fracture in Si</a></li>
<li class="toctree-l3"><a class="reference internal" href="adaptive-qmmm-step3.html">Step 3: LOTF hybrid MD simulation of fracture in Si</a></li>
<li class="toctree-l3 current"><a class="current reference internal" href="#">Solutions</a><ul>
<li class="toctree-l4"><a class="reference internal" href="#step-1-solution-make-crack-py">Step 1 solution — <code class="docutils literal notranslate"><span class="pre">make_crack.py</span></code></a></li>
<li class="toctree-l4"><a class="reference internal" href="#step-2-solution-run-crack-classical-py">Step 2 solution — <code class="docutils literal notranslate"><span class="pre">run_crack_classical.py</span></code></a></li>
<li class="toctree-l4"><a class="reference internal" href="#step-3-solution-run-crack-lotf-py">Step 3 solution — <code class="docutils literal notranslate"><span class="pre">run_crack_lotf.py</span></code></a></li>
</ul>
</li>
<li class="toctree-l3"><a class="reference internal" href="adaptive-qmmm-references.html">References</a></li>
</ul>
</li>
</ul>
</li>
<li class="toctree-l1"><a class="reference internal" href="../quippy.html">Quippy library reference</a></li>
<li class="toctree-l1"><a class="reference internal" href="../visualisation.html">Visualisation Tools</a></li>
<li class="toctree-l1"><a class="reference internal" href="../fortran_wrapper.html">Appendix: wrapping Fortran 90 Code</a></li>
<li class="toctree-l1"><a class="reference internal" href="../teach_sparse.html">Command line options of the teach_sparse main program</a></li>
</ul>

            
          
    <p></p>
    <ul>
        <li><a href="../genindex.html">Index</a></li>
        <li><a href="../py-modindex.html">Module Index</a></li>
        <li><a href="../search.html">Search Page</a></li>
        <li><a href="../_modules/index.html">Module Source Code</a></li>
    </ul>
  
        </div>
      </div>
    </nav>

    <section data-toggle="wy-nav-shift" class="wy-nav-content-wrap">

      
      <nav class="wy-nav-top" aria-label="top navigation">
        
          <i data-toggle="wy-nav-top" class="fa fa-bars"></i>
          <a href="../index.html">quippy</a>
        
      </nav>


      <div class="wy-nav-content">
        
        <div class="rst-content">
        
          















<div role="navigation" aria-label="breadcrumbs navigation">

  <ul class="wy-breadcrumbs">
    
      <li><a href="../index.html">Docs</a> &raquo;</li>
        
          <li><a href="index.html">Tutorials</a> &raquo;</li>
        
          <li><a href="adaptive-qmmm.html">Adaptive QM/MM MD of Fracture in Silicon</a> &raquo;</li>
        
      <li>Solutions</li>
    
    
      <li class="wy-breadcrumbs-aside">
        
            
            <a href="../_sources/Tutorials/adaptive-qmmm-solutions.rst.txt" rel="nofollow"> View page source</a>
          
        
      </li>
    
  </ul>

  
  <hr/>
</div>
          <div role="main" class="document" itemscope="itemscope" itemtype="http://schema.org/Article">
           <div itemprop="articleBody">
            
  
<style>
/* CSS overrides for sphinx_rtd_theme */

/* 24px margin */
.nbinput.nblast,
.nboutput.nblast {
    margin-bottom: 19px;  /* padding has already 5px */
}

/* ... except between code cells! */
.nblast + .nbinput {
    margin-top: -19px;
}

.admonition > p:before {
    margin-right: 4px;  /* make room for the exclamation icon */
}

/* Fix math alignment, see https://github.com/rtfd/sphinx_rtd_theme/pull/686 */
.math {
    text-align: unset;
}
</style>
<div class="section" id="solutions">
<h1>Solutions<a class="headerlink" href="#solutions" title="Permalink to this headline">¶</a></h1>
<div class="section" id="step-1-solution-make-crack-py">
<span id="make-crack"></span><h2>Step 1 solution — <code class="docutils literal notranslate"><span class="pre">make_crack.py</span></code><a class="headerlink" href="#step-1-solution-make-crack-py" title="Permalink to this headline">¶</a></h2>
<p>Download <a class="reference download internal" download="" href="../_downloads/df41e1fe21e67c537e56c9c0bac21586/make_crack.py"><code class="xref download docutils literal notranslate"><span class="pre">make_crack.py</span></code></a></p>
<div class="highlight-default notranslate"><div class="highlight"><pre><span></span><span class="sd">&quot;&quot;&quot;</span>
<span class="sd">make_crack.py</span>

<span class="sd">Script to generate a crack slab, and apply initial strain ramp</span>

<span class="sd">James Kermode &lt;james.kermode@kcl.ac.uk&gt;</span>
<span class="sd">January 2013</span>
<span class="sd">&quot;&quot;&quot;</span>
<span class="kn">from</span> <span class="nn">ase.structure</span> <span class="k">import</span> <span class="n">bulk</span>
<span class="kn">from</span> <span class="nn">ase.lattice.cubic</span> <span class="k">import</span> <span class="n">Diamond</span>
<span class="kn">from</span> <span class="nn">ase.constraints</span> <span class="k">import</span> <span class="n">FixAtoms</span>
<span class="kn">import</span> <span class="nn">ase.units</span> <span class="k">as</span> <span class="nn">units</span>

<span class="kn">from</span> <span class="nn">quippy</span> <span class="k">import</span> <span class="n">set_fortran_indexing</span>
<span class="kn">from</span> <span class="nn">quippy.atoms</span> <span class="k">import</span> <span class="n">Atoms</span>
<span class="kn">from</span> <span class="nn">quippy.potential</span> <span class="k">import</span> <span class="n">Potential</span><span class="p">,</span> <span class="n">Minim</span>
<span class="kn">from</span> <span class="nn">quippy.elasticity</span> <span class="k">import</span> <span class="n">youngs_modulus</span><span class="p">,</span> <span class="n">poisson_ratio</span>
<span class="kn">from</span> <span class="nn">quippy.io</span> <span class="k">import</span> <span class="n">write</span>
<span class="kn">from</span> <span class="nn">quippy.crack</span> <span class="k">import</span> <span class="p">(</span><span class="n">print_crack_system</span><span class="p">,</span>
                          <span class="n">G_to_strain</span><span class="p">,</span>
                          <span class="n">thin_strip_displacement_y</span><span class="p">,</span>
                          <span class="n">find_crack_tip_stress_field</span><span class="p">)</span>

<span class="c1"># ******* Start of parameters ***********</span>

<span class="c1"># There are three possible crack systems, choose one and uncomment it</span>

<span class="c1"># System 1. (111)[0-11]</span>
<span class="n">crack_direction</span> <span class="o">=</span> <span class="p">(</span><span class="o">-</span><span class="mi">2</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">)</span>      <span class="c1"># Miller index of x-axis</span>
<span class="n">cleavage_plane</span> <span class="o">=</span> <span class="p">(</span><span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">)</span>        <span class="c1"># Miller index of y-axis</span>
<span class="n">crack_front</span> <span class="o">=</span> <span class="p">(</span><span class="mi">0</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="o">-</span><span class="mi">1</span><span class="p">)</span>          <span class="c1"># Miller index of z-axis</span>

<span class="c1"># # System 2. (110)[001]</span>
<span class="c1"># crack_direction = (1,-1,0)</span>
<span class="c1"># cleavage_plane = (1,1,0)</span>
<span class="c1"># crack_front = (0,0,1)</span>

<span class="c1"># # System 3. (110)[1-10]</span>
<span class="c1"># crack_direction = (0,0,-1)</span>
<span class="c1"># cleavage_plane = (1,1,0)</span>
<span class="c1"># crack_front = (1,-1,0)</span>

<span class="n">width</span> <span class="o">=</span> <span class="mf">200.0</span><span class="o">*</span><span class="n">units</span><span class="o">.</span><span class="n">Ang</span>              <span class="c1"># Width of crack slab</span>
<span class="n">height</span> <span class="o">=</span> <span class="mf">100.0</span><span class="o">*</span><span class="n">units</span><span class="o">.</span><span class="n">Ang</span>             <span class="c1"># Height of crack slab</span>
<span class="n">vacuum</span> <span class="o">=</span> <span class="mf">100.0</span><span class="o">*</span><span class="n">units</span><span class="o">.</span><span class="n">Ang</span>             <span class="c1"># Amount of vacuum around slab</span>
<span class="n">crack_seed_length</span> <span class="o">=</span> <span class="mf">40.0</span><span class="o">*</span><span class="n">units</span><span class="o">.</span><span class="n">Ang</span>   <span class="c1"># Length of seed crack</span>
<span class="n">strain_ramp_length</span> <span class="o">=</span> <span class="mf">30.0</span><span class="o">*</span><span class="n">units</span><span class="o">.</span><span class="n">Ang</span>  <span class="c1"># Distance over which strain is ramped up</span>
<span class="n">initial_G</span> <span class="o">=</span> <span class="mf">5.0</span><span class="o">*</span><span class="p">(</span><span class="n">units</span><span class="o">.</span><span class="n">J</span><span class="o">/</span><span class="n">units</span><span class="o">.</span><span class="n">m</span><span class="o">**</span><span class="mi">2</span><span class="p">)</span> <span class="c1"># Initial energy flow to crack tip</span>

<span class="n">relax_fmax</span> <span class="o">=</span> <span class="mf">0.025</span><span class="o">*</span><span class="n">units</span><span class="o">.</span><span class="n">eV</span><span class="o">/</span><span class="n">units</span><span class="o">.</span><span class="n">Ang</span>  <span class="c1"># Maximum force criteria for relaxation</span>

<span class="n">param_file</span> <span class="o">=</span> <span class="s1">&#39;params.xml&#39;</span>            <span class="c1"># XML file containing interatomic potential parameters</span>
<span class="n">mm_init_args</span> <span class="o">=</span> <span class="s1">&#39;IP SW&#39;</span>               <span class="c1"># Initialisation arguments for the classical potential</span>

<span class="n">output_file</span> <span class="o">=</span> <span class="s1">&#39;crack.xyz&#39;</span>            <span class="c1"># File to which structure will be written</span>

<span class="c1"># ******* End of parameters *************</span>

<span class="n">set_fortran_indexing</span><span class="p">(</span><span class="kc">False</span><span class="p">)</span>

<span class="c1"># ********** Build unit cell ************</span>

<span class="c1"># 8-atom diamond cubic unit cell for silicon, with guess at lattice</span>
<span class="c1"># constant of 5.44 A</span>
<span class="n">si_bulk</span> <span class="o">=</span> <span class="n">bulk</span><span class="p">(</span><span class="s1">&#39;Si&#39;</span><span class="p">,</span> <span class="s1">&#39;diamond&#39;</span><span class="p">,</span> <span class="n">a</span><span class="o">=</span><span class="mf">5.44</span><span class="p">,</span> <span class="n">cubic</span><span class="o">=</span><span class="kc">True</span><span class="p">)</span>


<span class="c1"># ********** Setup potential ************</span>

<span class="c1"># Stillinger-Weber (SW) classical interatomic potential, from QUIP</span>
<span class="n">mm_pot</span> <span class="o">=</span> <span class="n">Potential</span><span class="p">(</span><span class="n">mm_init_args</span><span class="p">,</span>
                   <span class="n">param_filename</span><span class="o">=</span><span class="n">param_file</span><span class="p">)</span>


<span class="c1"># ***** Find eqm. lattice constant ******</span>

<span class="c1"># find the equilibrium lattice constant by minimising atoms wrt virial</span>
<span class="c1"># tensor given by SW pot (possibly replace this with parabola fit in</span>
<span class="c1"># another script and hardcoded a0 here)</span>
<span class="n">si_bulk</span><span class="o">.</span><span class="n">set_calculator</span><span class="p">(</span><span class="n">mm_pot</span><span class="p">)</span>

<span class="nb">print</span><span class="p">(</span><span class="s1">&#39;Minimising bulk unit cell...&#39;</span><span class="p">)</span>
<span class="n">minim</span> <span class="o">=</span> <span class="n">Minim</span><span class="p">(</span><span class="n">si_bulk</span><span class="p">,</span> <span class="n">relax_positions</span><span class="o">=</span><span class="kc">True</span><span class="p">,</span> <span class="n">relax_cell</span><span class="o">=</span><span class="kc">True</span><span class="p">)</span>
<span class="n">minim</span><span class="o">.</span><span class="n">run</span><span class="p">(</span><span class="n">fmax</span><span class="o">=</span><span class="mf">1e-4</span><span class="p">)</span>
    
<span class="n">a0</span> <span class="o">=</span> <span class="n">si_bulk</span><span class="o">.</span><span class="n">cell</span><span class="p">[</span><span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">]</span>
<span class="nb">print</span><span class="p">(</span><span class="s1">&#39;Lattice constant </span><span class="si">%.3f</span><span class="s1"> A</span><span class="se">\n</span><span class="s1">&#39;</span> <span class="o">%</span> <span class="n">a0</span><span class="p">)</span>

<span class="c1"># make a new bulk cell with correct a0 (so that off-diagonal lattice values are exactly zero)</span>
<span class="n">si_bulk</span> <span class="o">=</span> <span class="n">bulk</span><span class="p">(</span><span class="s1">&#39;Si&#39;</span><span class="p">,</span> <span class="s1">&#39;diamond&#39;</span><span class="p">,</span> <span class="n">a</span><span class="o">=</span><span class="n">a0</span><span class="p">,</span> <span class="n">cubic</span><span class="o">=</span><span class="kc">True</span><span class="p">)</span>
<span class="n">si_bulk</span><span class="o">.</span><span class="n">set_calculator</span><span class="p">(</span><span class="n">mm_pot</span><span class="p">)</span>

<span class="c1"># ******* Find elastic constants *******</span>

<span class="c1"># Get 6x6 matrix of elastic constants C_ij</span>
<span class="n">c</span> <span class="o">=</span> <span class="n">mm_pot</span><span class="o">.</span><span class="n">get_elastic_constants</span><span class="p">(</span><span class="n">si_bulk</span><span class="p">)</span>
<span class="nb">print</span><span class="p">(</span><span class="s1">&#39;Elastic constants (GPa):&#39;</span><span class="p">)</span>
<span class="nb">print</span><span class="p">((</span><span class="n">c</span> <span class="o">/</span> <span class="n">units</span><span class="o">.</span><span class="n">GPa</span><span class="p">)</span><span class="o">.</span><span class="n">round</span><span class="p">(</span><span class="mi">0</span><span class="p">))</span>
<span class="nb">print</span><span class="p">(</span><span class="s1">&#39;&#39;</span><span class="p">)</span>

<span class="n">E</span> <span class="o">=</span> <span class="n">youngs_modulus</span><span class="p">(</span><span class="n">c</span><span class="p">,</span> <span class="n">cleavage_plane</span><span class="p">)</span>
<span class="nb">print</span><span class="p">(</span><span class="s1">&#39;Young</span><span class="se">\&#39;</span><span class="s1">s modulus </span><span class="si">%.1f</span><span class="s1"> GPa&#39;</span> <span class="o">%</span> <span class="p">(</span><span class="n">E</span> <span class="o">/</span> <span class="n">units</span><span class="o">.</span><span class="n">GPa</span><span class="p">))</span>
<span class="n">nu</span> <span class="o">=</span> <span class="n">poisson_ratio</span><span class="p">(</span><span class="n">c</span><span class="p">,</span> <span class="n">cleavage_plane</span><span class="p">,</span> <span class="n">crack_direction</span><span class="p">)</span>
<span class="nb">print</span><span class="p">(</span><span class="s1">&#39;Poisson ratio </span><span class="si">% .3f</span><span class="se">\n</span><span class="s1">&#39;</span> <span class="o">%</span> <span class="n">nu</span><span class="p">)</span>

<span class="c1"># **** Setup crack slab unit cell ******</span>

<span class="n">print_crack_system</span><span class="p">(</span><span class="n">crack_direction</span><span class="p">,</span> <span class="n">cleavage_plane</span><span class="p">,</span> <span class="n">crack_front</span><span class="p">)</span>

<span class="c1"># now, we build system aligned with requested crystallographic orientation</span>
<span class="n">unit_slab</span> <span class="o">=</span> <span class="n">Diamond</span><span class="p">(</span><span class="n">directions</span><span class="o">=</span><span class="p">[</span><span class="n">crack_direction</span><span class="p">,</span>
                                <span class="n">cleavage_plane</span><span class="p">,</span>
                                <span class="n">crack_front</span><span class="p">],</span>
                    <span class="n">size</span><span class="o">=</span><span class="p">(</span><span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">),</span>
                    <span class="n">symbol</span><span class="o">=</span><span class="s1">&#39;Si&#39;</span><span class="p">,</span>
                    <span class="n">pbc</span><span class="o">=</span><span class="kc">True</span><span class="p">,</span>
                    <span class="n">latticeconstant</span><span class="o">=</span><span class="n">a0</span><span class="p">)</span>

<span class="c1"># You could check elastic constants of this rotated system:</span>
<span class="c1"># should lead to same Young&#39;s modulus and Poisson ratio</span>

<span class="nb">print</span><span class="p">(</span><span class="s1">&#39;Unit slab with </span><span class="si">%d</span><span class="s1"> atoms per unit cell:&#39;</span> <span class="o">%</span> <span class="nb">len</span><span class="p">(</span><span class="n">unit_slab</span><span class="p">))</span>
<span class="nb">print</span><span class="p">(</span><span class="n">unit_slab</span><span class="o">.</span><span class="n">cell</span><span class="p">)</span>
<span class="nb">print</span><span class="p">(</span><span class="s1">&#39;&#39;</span><span class="p">)</span>

<span class="c1"># center vertically half way along the vertical bond between atoms 0 and 1</span>
<span class="n">unit_slab</span><span class="o">.</span><span class="n">positions</span><span class="p">[:,</span> <span class="mi">1</span><span class="p">]</span> <span class="o">+=</span> <span class="p">(</span><span class="n">unit_slab</span><span class="o">.</span><span class="n">positions</span><span class="p">[</span><span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">]</span> <span class="o">-</span>
                              <span class="n">unit_slab</span><span class="o">.</span><span class="n">positions</span><span class="p">[</span><span class="mi">0</span><span class="p">,</span> <span class="mi">1</span><span class="p">])</span> <span class="o">/</span> <span class="mf">2.0</span>

<span class="c1"># map positions back into unit cell</span>
<span class="n">unit_slab</span><span class="o">.</span><span class="n">set_scaled_positions</span><span class="p">(</span><span class="n">unit_slab</span><span class="o">.</span><span class="n">get_scaled_positions</span><span class="p">())</span>

<span class="c1"># Make a surface unit cell by adding some vaccum along y</span>
<span class="n">surface</span> <span class="o">=</span> <span class="n">unit_slab</span><span class="o">.</span><span class="n">copy</span><span class="p">()</span>
<span class="n">surface</span><span class="o">.</span><span class="n">center</span><span class="p">(</span><span class="n">vacuum</span><span class="p">,</span> <span class="n">axis</span><span class="o">=</span><span class="mi">1</span><span class="p">)</span>


<span class="c1"># ********** Surface energy ************</span>

<span class="c1"># Calculate surface energy per unit area</span>
<span class="n">surface</span><span class="o">.</span><span class="n">set_calculator</span><span class="p">(</span><span class="n">mm_pot</span><span class="p">)</span>
<span class="n">E_surf</span> <span class="o">=</span> <span class="n">surface</span><span class="o">.</span><span class="n">get_potential_energy</span><span class="p">()</span>
<span class="n">E_per_atom_bulk</span> <span class="o">=</span> <span class="n">si_bulk</span><span class="o">.</span><span class="n">get_potential_energy</span><span class="p">()</span> <span class="o">/</span> <span class="nb">len</span><span class="p">(</span><span class="n">si_bulk</span><span class="p">)</span>
<span class="n">area</span> <span class="o">=</span> <span class="n">surface</span><span class="o">.</span><span class="n">get_volume</span><span class="p">()</span> <span class="o">/</span> <span class="n">surface</span><span class="o">.</span><span class="n">cell</span><span class="p">[</span><span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">]</span>
<span class="n">gamma</span> <span class="o">=</span> <span class="p">((</span><span class="n">E_surf</span> <span class="o">-</span> <span class="n">E_per_atom_bulk</span> <span class="o">*</span> <span class="nb">len</span><span class="p">(</span><span class="n">surface</span><span class="p">))</span> <span class="o">/</span>
         <span class="p">(</span><span class="mf">2.0</span> <span class="o">*</span> <span class="n">area</span><span class="p">))</span>

<span class="nb">print</span><span class="p">(</span><span class="s1">&#39;Surface energy of </span><span class="si">%s</span><span class="s1"> surface </span><span class="si">%.4f</span><span class="s1"> J/m^2</span><span class="se">\n</span><span class="s1">&#39;</span> <span class="o">%</span>
      <span class="p">(</span><span class="n">cleavage_plane</span><span class="p">,</span> <span class="n">gamma</span> <span class="o">/</span> <span class="p">(</span><span class="n">units</span><span class="o">.</span><span class="n">J</span> <span class="o">/</span> <span class="n">units</span><span class="o">.</span><span class="n">m</span> <span class="o">**</span> <span class="mi">2</span><span class="p">)))</span>


<span class="c1"># ***** Setup crack slab supercell *****</span>

<span class="c1"># Now we will build the full crack slab system,</span>
<span class="c1"># approximately matching requested width and height</span>
<span class="n">nx</span> <span class="o">=</span> <span class="nb">int</span><span class="p">(</span><span class="n">width</span> <span class="o">/</span> <span class="n">unit_slab</span><span class="o">.</span><span class="n">cell</span><span class="p">[</span><span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">])</span>
<span class="n">ny</span> <span class="o">=</span> <span class="nb">int</span><span class="p">(</span><span class="n">height</span> <span class="o">/</span> <span class="n">unit_slab</span><span class="o">.</span><span class="n">cell</span><span class="p">[</span><span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">])</span>

<span class="c1"># make sure ny is even so slab is centered on a bond</span>
<span class="k">if</span> <span class="n">ny</span> <span class="o">%</span> <span class="mi">2</span> <span class="o">==</span> <span class="mi">1</span><span class="p">:</span>
    <span class="n">ny</span> <span class="o">+=</span> <span class="mi">1</span>

<span class="c1"># make a supercell of unit_slab</span>
<span class="n">crack_slab</span> <span class="o">=</span> <span class="n">unit_slab</span> <span class="o">*</span> <span class="p">(</span><span class="n">nx</span><span class="p">,</span> <span class="n">ny</span><span class="p">,</span> <span class="mi">1</span><span class="p">)</span>

<span class="c1"># open up the cell along x and y by introducing some vaccum</span>
<span class="n">crack_slab</span><span class="o">.</span><span class="n">center</span><span class="p">(</span><span class="n">vacuum</span><span class="p">,</span> <span class="n">axis</span><span class="o">=</span><span class="mi">0</span><span class="p">)</span>
<span class="n">crack_slab</span><span class="o">.</span><span class="n">center</span><span class="p">(</span><span class="n">vacuum</span><span class="p">,</span> <span class="n">axis</span><span class="o">=</span><span class="mi">1</span><span class="p">)</span>

<span class="c1"># centre the slab on the origin</span>
<span class="n">crack_slab</span><span class="o">.</span><span class="n">positions</span><span class="p">[:,</span> <span class="mi">0</span><span class="p">]</span> <span class="o">-=</span> <span class="n">crack_slab</span><span class="o">.</span><span class="n">positions</span><span class="p">[:,</span> <span class="mi">0</span><span class="p">]</span><span class="o">.</span><span class="n">mean</span><span class="p">()</span>
<span class="n">crack_slab</span><span class="o">.</span><span class="n">positions</span><span class="p">[:,</span> <span class="mi">1</span><span class="p">]</span> <span class="o">-=</span> <span class="n">crack_slab</span><span class="o">.</span><span class="n">positions</span><span class="p">[:,</span> <span class="mi">1</span><span class="p">]</span><span class="o">.</span><span class="n">mean</span><span class="p">()</span>

<span class="n">orig_width</span> <span class="o">=</span> <span class="p">(</span><span class="n">crack_slab</span><span class="o">.</span><span class="n">positions</span><span class="p">[:,</span> <span class="mi">0</span><span class="p">]</span><span class="o">.</span><span class="n">max</span><span class="p">()</span> <span class="o">-</span>
              <span class="n">crack_slab</span><span class="o">.</span><span class="n">positions</span><span class="p">[:,</span> <span class="mi">0</span><span class="p">]</span><span class="o">.</span><span class="n">min</span><span class="p">())</span>
<span class="n">orig_height</span> <span class="o">=</span> <span class="p">(</span><span class="n">crack_slab</span><span class="o">.</span><span class="n">positions</span><span class="p">[:,</span> <span class="mi">1</span><span class="p">]</span><span class="o">.</span><span class="n">max</span><span class="p">()</span> <span class="o">-</span>
               <span class="n">crack_slab</span><span class="o">.</span><span class="n">positions</span><span class="p">[:,</span> <span class="mi">1</span><span class="p">]</span><span class="o">.</span><span class="n">min</span><span class="p">())</span>

<span class="nb">print</span><span class="p">((</span><span class="s1">&#39;Made slab with </span><span class="si">%d</span><span class="s1"> atoms, original width and height: </span><span class="si">%.1f</span><span class="s1"> x </span><span class="si">%.1f</span><span class="s1"> A^2&#39;</span> <span class="o">%</span>
       <span class="p">(</span><span class="nb">len</span><span class="p">(</span><span class="n">crack_slab</span><span class="p">),</span> <span class="n">orig_width</span><span class="p">,</span> <span class="n">orig_height</span><span class="p">)))</span>

<span class="n">top</span> <span class="o">=</span> <span class="n">crack_slab</span><span class="o">.</span><span class="n">positions</span><span class="p">[:,</span> <span class="mi">1</span><span class="p">]</span><span class="o">.</span><span class="n">max</span><span class="p">()</span>
<span class="n">bottom</span> <span class="o">=</span> <span class="n">crack_slab</span><span class="o">.</span><span class="n">positions</span><span class="p">[:,</span> <span class="mi">1</span><span class="p">]</span><span class="o">.</span><span class="n">min</span><span class="p">()</span>
<span class="n">left</span> <span class="o">=</span> <span class="n">crack_slab</span><span class="o">.</span><span class="n">positions</span><span class="p">[:,</span> <span class="mi">0</span><span class="p">]</span><span class="o">.</span><span class="n">min</span><span class="p">()</span>
<span class="n">right</span> <span class="o">=</span> <span class="n">crack_slab</span><span class="o">.</span><span class="n">positions</span><span class="p">[:,</span> <span class="mi">0</span><span class="p">]</span><span class="o">.</span><span class="n">max</span><span class="p">()</span>

<span class="c1"># fix atoms in the top and bottom rows</span>
<span class="n">fixed_mask</span> <span class="o">=</span> <span class="p">((</span><span class="nb">abs</span><span class="p">(</span><span class="n">crack_slab</span><span class="o">.</span><span class="n">positions</span><span class="p">[:,</span> <span class="mi">1</span><span class="p">]</span> <span class="o">-</span> <span class="n">top</span><span class="p">)</span> <span class="o">&lt;</span> <span class="mf">1.0</span><span class="p">)</span> <span class="o">|</span>
              <span class="p">(</span><span class="nb">abs</span><span class="p">(</span><span class="n">crack_slab</span><span class="o">.</span><span class="n">positions</span><span class="p">[:,</span> <span class="mi">1</span><span class="p">]</span> <span class="o">-</span> <span class="n">bottom</span><span class="p">)</span> <span class="o">&lt;</span> <span class="mf">1.0</span><span class="p">))</span>
<span class="n">const</span> <span class="o">=</span> <span class="n">FixAtoms</span><span class="p">(</span><span class="n">mask</span><span class="o">=</span><span class="n">fixed_mask</span><span class="p">)</span>
<span class="n">crack_slab</span><span class="o">.</span><span class="n">set_constraint</span><span class="p">(</span><span class="n">const</span><span class="p">)</span>
<span class="nb">print</span><span class="p">(</span><span class="s1">&#39;Fixed </span><span class="si">%d</span><span class="s1"> atoms</span><span class="se">\n</span><span class="s1">&#39;</span> <span class="o">%</span> <span class="n">fixed_mask</span><span class="o">.</span><span class="n">sum</span><span class="p">())</span>


<span class="c1"># ****** Apply initial strain ramp *****</span>

<span class="n">strain</span> <span class="o">=</span> <span class="n">G_to_strain</span><span class="p">(</span><span class="n">initial_G</span><span class="p">,</span> <span class="n">E</span><span class="p">,</span> <span class="n">nu</span><span class="p">,</span> <span class="n">orig_height</span><span class="p">)</span>

<span class="n">crack_slab</span><span class="o">.</span><span class="n">positions</span><span class="p">[:,</span> <span class="mi">1</span><span class="p">]</span> <span class="o">+=</span> <span class="n">thin_strip_displacement_y</span><span class="p">(</span>
                                 <span class="n">crack_slab</span><span class="o">.</span><span class="n">positions</span><span class="p">[:,</span> <span class="mi">0</span><span class="p">],</span>
                                 <span class="n">crack_slab</span><span class="o">.</span><span class="n">positions</span><span class="p">[:,</span> <span class="mi">1</span><span class="p">],</span>
                                 <span class="n">strain</span><span class="p">,</span>
                                 <span class="n">left</span> <span class="o">+</span> <span class="n">crack_seed_length</span><span class="p">,</span>
                                 <span class="n">left</span> <span class="o">+</span> <span class="n">crack_seed_length</span> <span class="o">+</span> <span class="n">strain_ramp_length</span><span class="p">)</span>

<span class="nb">print</span><span class="p">(</span><span class="s1">&#39;Applied initial load: strain=</span><span class="si">%.4f</span><span class="s1">, G=</span><span class="si">%.2f</span><span class="s1"> J/m^2&#39;</span> <span class="o">%</span>
      <span class="p">(</span><span class="n">strain</span><span class="p">,</span> <span class="n">initial_G</span> <span class="o">/</span> <span class="p">(</span><span class="n">units</span><span class="o">.</span><span class="n">J</span> <span class="o">/</span> <span class="n">units</span><span class="o">.</span><span class="n">m</span><span class="o">**</span><span class="mi">2</span><span class="p">)))</span>


<span class="c1"># ***** Relaxation of crack slab  *****</span>

<span class="c1"># optionally, relax the slab, keeping top and bottom rows fixed</span>
<span class="nb">print</span><span class="p">(</span><span class="s1">&#39;Relaxing slab...&#39;</span><span class="p">)</span>
<span class="n">crack_slab</span><span class="o">.</span><span class="n">set_calculator</span><span class="p">(</span><span class="n">mm_pot</span><span class="p">)</span>
<span class="n">minim</span> <span class="o">=</span> <span class="n">Minim</span><span class="p">(</span><span class="n">crack_slab</span><span class="p">,</span> <span class="n">relax_positions</span><span class="o">=</span><span class="kc">True</span><span class="p">,</span> <span class="n">relax_cell</span><span class="o">=</span><span class="kc">False</span><span class="p">)</span>
<span class="n">minim</span><span class="o">.</span><span class="n">run</span><span class="p">(</span><span class="n">fmax</span><span class="o">=</span><span class="n">relax_fmax</span><span class="p">)</span>

<span class="c1"># Find initial position of crack tip</span>
<span class="n">crack_pos</span> <span class="o">=</span> <span class="n">find_crack_tip_stress_field</span><span class="p">(</span><span class="n">crack_slab</span><span class="p">,</span> <span class="n">calc</span><span class="o">=</span><span class="n">mm_pot</span><span class="p">)</span>
<span class="nb">print</span> <span class="s1">&#39;Found crack tip at position </span><span class="si">%s</span><span class="s1">&#39;</span> <span class="o">%</span> <span class="n">crack_pos</span>

<span class="c1"># Save all calculated materials properties inside the Atoms object</span>
<span class="n">crack_slab</span><span class="o">.</span><span class="n">info</span><span class="p">[</span><span class="s1">&#39;nneightol&#39;</span><span class="p">]</span> <span class="o">=</span> <span class="mf">1.3</span> <span class="c1"># nearest neighbour tolerance</span>
<span class="n">crack_slab</span><span class="o">.</span><span class="n">info</span><span class="p">[</span><span class="s1">&#39;LatticeConstant&#39;</span><span class="p">]</span> <span class="o">=</span> <span class="n">a0</span>
<span class="n">crack_slab</span><span class="o">.</span><span class="n">info</span><span class="p">[</span><span class="s1">&#39;C11&#39;</span><span class="p">]</span> <span class="o">=</span> <span class="n">c</span><span class="p">[</span><span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">]</span>
<span class="n">crack_slab</span><span class="o">.</span><span class="n">info</span><span class="p">[</span><span class="s1">&#39;C12&#39;</span><span class="p">]</span> <span class="o">=</span> <span class="n">c</span><span class="p">[</span><span class="mi">0</span><span class="p">,</span> <span class="mi">1</span><span class="p">]</span>
<span class="n">crack_slab</span><span class="o">.</span><span class="n">info</span><span class="p">[</span><span class="s1">&#39;C44&#39;</span><span class="p">]</span> <span class="o">=</span> <span class="n">c</span><span class="p">[</span><span class="mi">3</span><span class="p">,</span> <span class="mi">3</span><span class="p">]</span>
<span class="n">crack_slab</span><span class="o">.</span><span class="n">info</span><span class="p">[</span><span class="s1">&#39;YoungsModulus&#39;</span><span class="p">]</span> <span class="o">=</span> <span class="n">E</span>
<span class="n">crack_slab</span><span class="o">.</span><span class="n">info</span><span class="p">[</span><span class="s1">&#39;PoissonRatio_yx&#39;</span><span class="p">]</span> <span class="o">=</span> <span class="n">nu</span>
<span class="n">crack_slab</span><span class="o">.</span><span class="n">info</span><span class="p">[</span><span class="s1">&#39;SurfaceEnergy&#39;</span><span class="p">]</span> <span class="o">=</span> <span class="n">gamma</span>
<span class="n">crack_slab</span><span class="o">.</span><span class="n">info</span><span class="p">[</span><span class="s1">&#39;OrigWidth&#39;</span><span class="p">]</span> <span class="o">=</span> <span class="n">orig_width</span>
<span class="n">crack_slab</span><span class="o">.</span><span class="n">info</span><span class="p">[</span><span class="s1">&#39;OrigHeight&#39;</span><span class="p">]</span> <span class="o">=</span> <span class="n">orig_height</span>
<span class="n">crack_slab</span><span class="o">.</span><span class="n">info</span><span class="p">[</span><span class="s1">&#39;CrackDirection&#39;</span><span class="p">]</span> <span class="o">=</span> <span class="n">crack_direction</span>
<span class="n">crack_slab</span><span class="o">.</span><span class="n">info</span><span class="p">[</span><span class="s1">&#39;CleavagePlane&#39;</span><span class="p">]</span> <span class="o">=</span> <span class="n">cleavage_plane</span>
<span class="n">crack_slab</span><span class="o">.</span><span class="n">info</span><span class="p">[</span><span class="s1">&#39;CrackFront&#39;</span><span class="p">]</span> <span class="o">=</span> <span class="n">crack_front</span>
<span class="n">crack_slab</span><span class="o">.</span><span class="n">info</span><span class="p">[</span><span class="s1">&#39;strain&#39;</span><span class="p">]</span> <span class="o">=</span> <span class="n">strain</span>
<span class="n">crack_slab</span><span class="o">.</span><span class="n">info</span><span class="p">[</span><span class="s1">&#39;G&#39;</span><span class="p">]</span> <span class="o">=</span> <span class="n">initial_G</span>
<span class="n">crack_slab</span><span class="o">.</span><span class="n">info</span><span class="p">[</span><span class="s1">&#39;CrackPos&#39;</span><span class="p">]</span> <span class="o">=</span> <span class="n">crack_pos</span>
<span class="n">crack_slab</span><span class="o">.</span><span class="n">info</span><span class="p">[</span><span class="s1">&#39;is_cracked&#39;</span><span class="p">]</span> <span class="o">=</span> <span class="kc">False</span>


<span class="c1"># ******** Save output file **********</span>

<span class="c1"># save results in extended XYZ format, including extra properties and info</span>
<span class="nb">print</span><span class="p">(</span><span class="s1">&#39;Writing crack slab to file </span><span class="si">%s</span><span class="s1">&#39;</span> <span class="o">%</span> <span class="n">output_file</span><span class="p">)</span>
<span class="n">write</span><span class="p">(</span><span class="n">output_file</span><span class="p">,</span> <span class="n">crack_slab</span><span class="p">)</span>
</pre></div>
</div>
</div>
<div class="section" id="step-2-solution-run-crack-classical-py">
<span id="run-crack-classical"></span><h2>Step 2 solution — <code class="docutils literal notranslate"><span class="pre">run_crack_classical.py</span></code><a class="headerlink" href="#step-2-solution-run-crack-classical-py" title="Permalink to this headline">¶</a></h2>
<p>Download <a class="reference download internal" download="" href="../_downloads/97a9026d6ccedee9572afdff54b93c12/run_crack_classical.py"><code class="xref download docutils literal notranslate"><span class="pre">run_crack_classical.py</span></code></a></p>
<div class="highlight-default notranslate"><div class="highlight"><pre><span></span><span class="sd">&quot;&quot;&quot;</span>
<span class="sd">run_crack_classical.py</span>

<span class="sd">Script to run classical molecular dynamics for a crack slab,</span>
<span class="sd">incrementing the load in small steps until fracture starts.</span>

<span class="sd">James Kermode &lt;james.kermode@kcl.ac.uk&gt;</span>
<span class="sd">January 2013</span>
<span class="sd">&quot;&quot;&quot;</span>

<span class="kn">import</span> <span class="nn">numpy</span> <span class="k">as</span> <span class="nn">np</span>

<span class="kn">from</span> <span class="nn">ase.constraints</span> <span class="k">import</span> <span class="n">FixAtoms</span>
<span class="kn">from</span> <span class="nn">ase.md.verlet</span> <span class="k">import</span> <span class="n">VelocityVerlet</span>
<span class="kn">from</span> <span class="nn">ase.md.velocitydistribution</span> <span class="k">import</span> <span class="n">MaxwellBoltzmannDistribution</span>
<span class="kn">import</span> <span class="nn">ase.units</span> <span class="k">as</span> <span class="nn">units</span>

<span class="kn">from</span> <span class="nn">quippy</span> <span class="k">import</span> <span class="n">set_fortran_indexing</span>
<span class="kn">from</span> <span class="nn">quippy.atoms</span> <span class="k">import</span> <span class="n">Atoms</span>
<span class="kn">from</span> <span class="nn">quippy.potential</span> <span class="k">import</span> <span class="n">Potential</span>
<span class="kn">from</span> <span class="nn">quippy.io</span> <span class="k">import</span> <span class="n">AtomsWriter</span>

<span class="kn">from</span> <span class="nn">quippy.crack</span> <span class="k">import</span> <span class="p">(</span><span class="n">get_strain</span><span class="p">,</span>
                          <span class="n">get_energy_release_rate</span><span class="p">,</span>
                          <span class="n">ConstantStrainRate</span><span class="p">,</span>
                          <span class="n">find_crack_tip_stress_field</span><span class="p">)</span>

<span class="c1"># ******* Start of parameters ***********</span>

<span class="n">input_file</span> <span class="o">=</span> <span class="s1">&#39;crack.xyz&#39;</span>         <span class="c1"># File from which to read crack slab structure</span>
<span class="n">sim_T</span> <span class="o">=</span> <span class="mf">300.0</span><span class="o">*</span><span class="n">units</span><span class="o">.</span><span class="n">kB</span>           <span class="c1"># Simulation temperature</span>
<span class="n">nsteps</span> <span class="o">=</span> <span class="mi">10000</span>                   <span class="c1"># Total number of timesteps to run for</span>
<span class="n">timestep</span> <span class="o">=</span> <span class="mf">1.0</span><span class="o">*</span><span class="n">units</span><span class="o">.</span><span class="n">fs</span>          <span class="c1"># Timestep (NB: time base units are not fs!)</span>
<span class="n">cutoff_skin</span> <span class="o">=</span> <span class="mf">2.0</span><span class="o">*</span><span class="n">units</span><span class="o">.</span><span class="n">Ang</span>      <span class="c1"># Amount by which potential cutoff is increased</span>
                                 <span class="c1"># for neighbour calculations</span>
<span class="n">tip_move_tol</span> <span class="o">=</span> <span class="mf">10.0</span>              <span class="c1"># Distance tip has to move before crack</span>
                                 <span class="c1"># is taken to be running</span>
<span class="n">strain_rate</span> <span class="o">=</span> <span class="mf">1e-5</span><span class="o">*</span><span class="p">(</span><span class="mi">1</span><span class="o">/</span><span class="n">units</span><span class="o">.</span><span class="n">fs</span><span class="p">)</span>  <span class="c1"># Strain rate</span>
<span class="n">traj_file</span> <span class="o">=</span> <span class="s1">&#39;traj.xyz&#39;</span>           <span class="c1"># Trajectory output file (NetCDF format)</span>
<span class="n">traj_interval</span> <span class="o">=</span> <span class="mi">10</span>               <span class="c1"># Number of time steps between</span>
                                 <span class="c1"># writing output frames</span>
<span class="n">param_file</span> <span class="o">=</span> <span class="s1">&#39;params.xml&#39;</span>        <span class="c1"># Filename of XML file containing</span>
                                 <span class="c1"># potential parameters</span>
<span class="n">mm_init_args</span> <span class="o">=</span> <span class="s1">&#39;IP SW&#39;</span>           <span class="c1"># Initialisation arguments for</span>
                                 <span class="c1"># classical potential</span>

<span class="c1"># ******* End of parameters *************</span>

<span class="n">set_fortran_indexing</span><span class="p">(</span><span class="kc">False</span><span class="p">)</span>

<span class="c1"># ********** Read input file ************</span>

<span class="nb">print</span> <span class="s1">&#39;Loading atoms from file </span><span class="si">%s</span><span class="s1">&#39;</span> <span class="o">%</span> <span class="n">input_file</span>
<span class="n">atoms</span> <span class="o">=</span> <span class="n">Atoms</span><span class="p">(</span><span class="n">input_file</span><span class="p">)</span>

<span class="n">orig_height</span> <span class="o">=</span> <span class="n">atoms</span><span class="o">.</span><span class="n">info</span><span class="p">[</span><span class="s1">&#39;OrigHeight&#39;</span><span class="p">]</span>
<span class="n">orig_crack_pos</span> <span class="o">=</span> <span class="n">atoms</span><span class="o">.</span><span class="n">info</span><span class="p">[</span><span class="s1">&#39;CrackPos&#39;</span><span class="p">]</span><span class="o">.</span><span class="n">copy</span><span class="p">()</span>

<span class="c1"># ***** Setup constraints *******</span>

<span class="n">top</span> <span class="o">=</span> <span class="n">atoms</span><span class="o">.</span><span class="n">positions</span><span class="p">[:,</span> <span class="mi">1</span><span class="p">]</span><span class="o">.</span><span class="n">max</span><span class="p">()</span>
<span class="n">bottom</span> <span class="o">=</span> <span class="n">atoms</span><span class="o">.</span><span class="n">positions</span><span class="p">[:,</span> <span class="mi">1</span><span class="p">]</span><span class="o">.</span><span class="n">min</span><span class="p">()</span>
<span class="n">left</span> <span class="o">=</span> <span class="n">atoms</span><span class="o">.</span><span class="n">positions</span><span class="p">[:,</span> <span class="mi">0</span><span class="p">]</span><span class="o">.</span><span class="n">min</span><span class="p">()</span>
<span class="n">right</span> <span class="o">=</span> <span class="n">atoms</span><span class="o">.</span><span class="n">positions</span><span class="p">[:,</span> <span class="mi">0</span><span class="p">]</span><span class="o">.</span><span class="n">max</span><span class="p">()</span>

<span class="c1"># fix atoms in the top and bottom rows</span>
<span class="n">fixed_mask</span> <span class="o">=</span> <span class="p">((</span><span class="nb">abs</span><span class="p">(</span><span class="n">atoms</span><span class="o">.</span><span class="n">positions</span><span class="p">[:,</span> <span class="mi">1</span><span class="p">]</span> <span class="o">-</span> <span class="n">top</span><span class="p">)</span> <span class="o">&lt;</span> <span class="mf">1.0</span><span class="p">)</span> <span class="o">|</span>
              <span class="p">(</span><span class="nb">abs</span><span class="p">(</span><span class="n">atoms</span><span class="o">.</span><span class="n">positions</span><span class="p">[:,</span> <span class="mi">1</span><span class="p">]</span> <span class="o">-</span> <span class="n">bottom</span><span class="p">)</span> <span class="o">&lt;</span> <span class="mf">1.0</span><span class="p">))</span>
<span class="n">fix_atoms</span> <span class="o">=</span> <span class="n">FixAtoms</span><span class="p">(</span><span class="n">mask</span><span class="o">=</span><span class="n">fixed_mask</span><span class="p">)</span>
<span class="nb">print</span><span class="p">(</span><span class="s1">&#39;Fixed </span><span class="si">%d</span><span class="s1"> atoms</span><span class="se">\n</span><span class="s1">&#39;</span> <span class="o">%</span> <span class="n">fixed_mask</span><span class="o">.</span><span class="n">sum</span><span class="p">())</span>
<span class="n">atoms</span><span class="o">.</span><span class="n">set_constraint</span><span class="p">([</span><span class="n">fix_atoms</span><span class="p">])</span>

<span class="c1"># Increase epsilon_yy applied to all atoms at constant strain rate</span>
<span class="n">strain_atoms</span> <span class="o">=</span> <span class="n">ConstantStrainRate</span><span class="p">(</span><span class="n">orig_height</span><span class="p">,</span> <span class="n">strain_rate</span><span class="o">*</span><span class="n">timestep</span><span class="p">)</span>

<span class="c1"># ******* Set up potentials and calculators ********</span>

<span class="n">mm_pot</span> <span class="o">=</span> <span class="n">Potential</span><span class="p">(</span><span class="n">mm_init_args</span><span class="p">,</span>
                   <span class="n">param_filename</span><span class="o">=</span><span class="n">param_file</span><span class="p">,</span>
                   <span class="n">cutoff_skin</span><span class="o">=</span><span class="n">cutoff_skin</span><span class="p">)</span>

<span class="n">atoms</span><span class="o">.</span><span class="n">set_calculator</span><span class="p">(</span><span class="n">mm_pot</span><span class="p">)</span>

<span class="c1"># ********* Setup and run MD ***********</span>

<span class="c1"># Set the initial temperature to 2*simT: it will then equilibriate to</span>
<span class="c1"># simT, by the virial theorem</span>
<span class="n">MaxwellBoltzmannDistribution</span><span class="p">(</span><span class="n">atoms</span><span class="p">,</span> <span class="mf">2.0</span><span class="o">*</span><span class="n">sim_T</span><span class="p">)</span>

<span class="c1"># Initialise the dynamical system</span>
<span class="n">dynamics</span> <span class="o">=</span> <span class="n">VelocityVerlet</span><span class="p">(</span><span class="n">atoms</span><span class="p">,</span> <span class="n">timestep</span><span class="p">)</span>

<span class="c1"># Print some information every time step</span>
<span class="k">def</span> <span class="nf">printstatus</span><span class="p">():</span>
    <span class="k">if</span> <span class="n">dynamics</span><span class="o">.</span><span class="n">nsteps</span> <span class="o">==</span> <span class="mi">1</span><span class="p">:</span>
        <span class="nb">print</span> <span class="s2">&quot;&quot;&quot;</span>
<span class="s2">State      Time/fs    Temp/K     Strain      G/(J/m^2)  CrackPos/A D(CrackPos)/A</span>
<span class="s2">---------------------------------------------------------------------------------&quot;&quot;&quot;</span>

    <span class="n">log_format</span> <span class="o">=</span> <span class="p">(</span><span class="s1">&#39;</span><span class="si">%(label)-4s%(time)12.1f%(temperature)12.6f</span><span class="s1">&#39;</span><span class="o">+</span>
                  <span class="s1">&#39;</span><span class="si">%(strain)12.5f%(G)12.4f%(crack_pos_x)12.2f</span><span class="s1">    (</span><span class="si">%(d_crack_pos_x)+5.2f</span><span class="s1">)&#39;</span><span class="p">)</span>

    <span class="n">atoms</span><span class="o">.</span><span class="n">info</span><span class="p">[</span><span class="s1">&#39;label&#39;</span><span class="p">]</span> <span class="o">=</span> <span class="s1">&#39;D&#39;</span>                  <span class="c1"># Label for the status line</span>
    <span class="n">atoms</span><span class="o">.</span><span class="n">info</span><span class="p">[</span><span class="s1">&#39;time&#39;</span><span class="p">]</span> <span class="o">=</span> <span class="n">dynamics</span><span class="o">.</span><span class="n">get_time</span><span class="p">()</span><span class="o">/</span><span class="n">units</span><span class="o">.</span><span class="n">fs</span>
    <span class="n">atoms</span><span class="o">.</span><span class="n">info</span><span class="p">[</span><span class="s1">&#39;temperature&#39;</span><span class="p">]</span> <span class="o">=</span> <span class="p">(</span><span class="n">atoms</span><span class="o">.</span><span class="n">get_kinetic_energy</span><span class="p">()</span> <span class="o">/</span>
                                 <span class="p">(</span><span class="mf">1.5</span><span class="o">*</span><span class="n">units</span><span class="o">.</span><span class="n">kB</span><span class="o">*</span><span class="nb">len</span><span class="p">(</span><span class="n">atoms</span><span class="p">)))</span>
    <span class="n">atoms</span><span class="o">.</span><span class="n">info</span><span class="p">[</span><span class="s1">&#39;strain&#39;</span><span class="p">]</span> <span class="o">=</span> <span class="n">get_strain</span><span class="p">(</span><span class="n">atoms</span><span class="p">)</span>
    <span class="n">atoms</span><span class="o">.</span><span class="n">info</span><span class="p">[</span><span class="s1">&#39;G&#39;</span><span class="p">]</span> <span class="o">=</span> <span class="n">get_energy_release_rate</span><span class="p">(</span><span class="n">atoms</span><span class="p">)</span><span class="o">/</span><span class="p">(</span><span class="n">units</span><span class="o">.</span><span class="n">J</span><span class="o">/</span><span class="n">units</span><span class="o">.</span><span class="n">m</span><span class="o">**</span><span class="mi">2</span><span class="p">)</span>

    <span class="n">crack_pos</span> <span class="o">=</span> <span class="n">find_crack_tip_stress_field</span><span class="p">(</span><span class="n">atoms</span><span class="p">,</span> <span class="n">calc</span><span class="o">=</span><span class="n">mm_pot</span><span class="p">)</span>
    <span class="n">atoms</span><span class="o">.</span><span class="n">info</span><span class="p">[</span><span class="s1">&#39;crack_pos_x&#39;</span><span class="p">]</span> <span class="o">=</span> <span class="n">crack_pos</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span>
    <span class="n">atoms</span><span class="o">.</span><span class="n">info</span><span class="p">[</span><span class="s1">&#39;d_crack_pos_x&#39;</span><span class="p">]</span> <span class="o">=</span> <span class="n">crack_pos</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span> <span class="o">-</span> <span class="n">orig_crack_pos</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span>

    <span class="nb">print</span> <span class="n">log_format</span> <span class="o">%</span> <span class="n">atoms</span><span class="o">.</span><span class="n">info</span>


<span class="n">dynamics</span><span class="o">.</span><span class="n">attach</span><span class="p">(</span><span class="n">printstatus</span><span class="p">)</span>

<span class="k">def</span> <span class="nf">atom_straining</span><span class="p">(</span><span class="n">atoms</span><span class="p">):</span>
    <span class="n">crack_pos</span> <span class="o">=</span> <span class="n">find_crack_tip_stress_field</span><span class="p">(</span><span class="n">atoms</span><span class="p">,</span> <span class="n">calc</span><span class="o">=</span><span class="n">mm_pot</span><span class="p">)</span>
    <span class="c1"># keep straining until the crack tip has advanced to tip_move_tol</span>
    <span class="k">if</span> <span class="ow">not</span> <span class="n">atoms</span><span class="o">.</span><span class="n">info</span><span class="p">[</span><span class="s1">&#39;is_cracked&#39;</span><span class="p">]</span> <span class="ow">and</span> <span class="p">(</span><span class="n">crack_pos</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span> <span class="o">-</span> <span class="n">orig_crack_pos</span><span class="p">[</span><span class="mi">0</span><span class="p">])</span> <span class="o">&lt;</span> <span class="n">tip_move_tol</span><span class="p">:</span>
      <span class="n">strain_atoms</span><span class="o">.</span><span class="n">apply_strain</span><span class="p">(</span><span class="n">atoms</span><span class="p">)</span>
    <span class="k">elif</span> <span class="ow">not</span> <span class="n">atoms</span><span class="o">.</span><span class="n">info</span><span class="p">[</span><span class="s1">&#39;is_cracked&#39;</span><span class="p">]:</span>
      <span class="n">atoms</span><span class="o">.</span><span class="n">info</span><span class="p">[</span><span class="s1">&#39;is_cracked&#39;</span><span class="p">]</span> <span class="o">=</span> <span class="kc">True</span>

<span class="n">dynamics</span><span class="o">.</span><span class="n">attach</span><span class="p">(</span><span class="n">atom_straining</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="n">atoms</span><span class="p">)</span>

<span class="c1"># Save frames to the trajectory every `traj_interval` time steps</span>
<span class="n">trajectory</span> <span class="o">=</span> <span class="n">AtomsWriter</span><span class="p">(</span><span class="n">traj_file</span><span class="p">)</span>
<span class="n">dynamics</span><span class="o">.</span><span class="n">attach</span><span class="p">(</span><span class="n">trajectory</span><span class="p">,</span> <span class="n">traj_interval</span><span class="p">,</span> <span class="n">atoms</span><span class="p">)</span>

<span class="c1"># Start running!</span>
<span class="n">dynamics</span><span class="o">.</span><span class="n">run</span><span class="p">(</span><span class="n">nsteps</span><span class="p">)</span>
</pre></div>
</div>
</div>
<div class="section" id="step-3-solution-run-crack-lotf-py">
<span id="run-crack-lotf"></span><h2>Step 3 solution — <code class="docutils literal notranslate"><span class="pre">run_crack_lotf.py</span></code><a class="headerlink" href="#step-3-solution-run-crack-lotf-py" title="Permalink to this headline">¶</a></h2>
<p>Download <a class="reference download internal" download="" href="../_downloads/6796d4b41d5c4d2b74d531df497da8ef/run_crack_lotf.py"><code class="xref download docutils literal notranslate"><span class="pre">run_crack_lotf.py</span></code></a></p>
<div class="highlight-default notranslate"><div class="highlight"><pre><span></span><span class="sd">&quot;&quot;&quot;</span>
<span class="sd">run_crack_lotf.py</span>
<span class="sd">Script to run LOTF molecular dynamics for a crack slab,</span>
<span class="sd">incrementing the load in small steps until fracture starts.</span>
<span class="sd">James Kermode &lt;james.kermode@kcl.ac.uk&gt;</span>
<span class="sd">February 2013</span>
<span class="sd">&quot;&quot;&quot;</span>
<span class="kn">import</span> <span class="nn">numpy</span> <span class="k">as</span> <span class="nn">np</span>

<span class="kn">from</span> <span class="nn">ase.constraints</span> <span class="k">import</span> <span class="n">FixAtoms</span>
<span class="kn">from</span> <span class="nn">ase.md.verlet</span> <span class="k">import</span> <span class="n">VelocityVerlet</span>
<span class="kn">from</span> <span class="nn">ase.md.velocitydistribution</span> <span class="k">import</span> <span class="n">MaxwellBoltzmannDistribution</span>
<span class="kn">import</span> <span class="nn">ase.units</span> <span class="k">as</span> <span class="nn">units</span>

<span class="kn">from</span> <span class="nn">quippy</span> <span class="k">import</span> <span class="n">set_fortran_indexing</span>
<span class="kn">from</span> <span class="nn">quippy.atoms</span> <span class="k">import</span> <span class="n">Atoms</span>
<span class="kn">from</span> <span class="nn">quippy.potential</span> <span class="k">import</span> <span class="n">Potential</span>
<span class="kn">from</span> <span class="nn">quippy.io</span> <span class="k">import</span> <span class="n">AtomsWriter</span>

<span class="kn">from</span> <span class="nn">quippy.crack</span> <span class="k">import</span> <span class="p">(</span><span class="n">get_strain</span><span class="p">,</span>
                          <span class="n">get_energy_release_rate</span><span class="p">,</span>
                          <span class="n">ConstantStrainRate</span><span class="p">,</span>
                          <span class="n">find_crack_tip_stress_field</span><span class="p">)</span>

<span class="c1"># additional requirements for the QM/MM simulation:</span>
<span class="kn">from</span> <span class="nn">quippy.potential</span> <span class="k">import</span> <span class="n">ForceMixingPotential</span>
<span class="kn">from</span> <span class="nn">quippy.lotf</span> <span class="k">import</span> <span class="n">LOTFDynamics</span><span class="p">,</span> <span class="n">update_hysteretic_qm_region</span>


<span class="c1"># ******* Start of parameters ***********</span>

<span class="n">input_file</span> <span class="o">=</span> <span class="s1">&#39;crack.xyz&#39;</span>         <span class="c1"># File from which to read crack slab structure</span>
<span class="n">sim_T</span> <span class="o">=</span> <span class="mf">300.0</span><span class="o">*</span><span class="n">units</span><span class="o">.</span><span class="n">kB</span>           <span class="c1"># Simulation temperature</span>
<span class="n">nsteps</span> <span class="o">=</span> <span class="mi">1000000</span>                   <span class="c1"># Total number of timesteps to run for</span>
<span class="n">timestep</span> <span class="o">=</span> <span class="mf">1.0</span><span class="o">*</span><span class="n">units</span><span class="o">.</span><span class="n">fs</span>          <span class="c1"># Timestep (NB: time base units are not fs!)</span>
<span class="n">cutoff_skin</span> <span class="o">=</span> <span class="mf">2.0</span><span class="o">*</span><span class="n">units</span><span class="o">.</span><span class="n">Ang</span>      <span class="c1"># Amount by which potential cutoff is increased</span>
                                 <span class="c1"># for neighbour calculations</span>
<span class="n">tip_move_tol</span> <span class="o">=</span> <span class="mf">10.0</span>              <span class="c1"># Distance tip has to move before crack</span>
                                 <span class="c1"># is taken to be running</span>
<span class="n">strain_rate</span> <span class="o">=</span> <span class="mf">1e-5</span><span class="o">*</span><span class="p">(</span><span class="mi">1</span><span class="o">/</span><span class="n">units</span><span class="o">.</span><span class="n">fs</span><span class="p">)</span>  <span class="c1"># Strain rate</span>
<span class="n">traj_file</span> <span class="o">=</span> <span class="s1">&#39;traj.xyz&#39;</span>           <span class="c1"># Trajectory output file (NetCDF or XYZ format)</span>
<span class="n">traj_interval</span> <span class="o">=</span> <span class="mi">10</span>               <span class="c1"># Number of time steps between</span>
                                 <span class="c1"># writing output frames</span>
<span class="n">param_file</span> <span class="o">=</span> <span class="s1">&#39;params.xml&#39;</span>        <span class="c1"># Filename of XML file containing</span>
                                 <span class="c1"># potential parameters</span>
<span class="n">mm_init_args</span> <span class="o">=</span> <span class="s1">&#39;IP SW&#39;</span>           <span class="c1"># Initialisation arguments for</span>
                                 <span class="c1"># classical potential</span>

<span class="c1"># additional parameters for the QM/MM simulation:</span>
<span class="n">qm_init_args</span> <span class="o">=</span> <span class="s1">&#39;TB DFTB&#39;</span>         <span class="c1"># Initialisation arguments for QM potential</span>
<span class="n">qm_inner_radius</span> <span class="o">=</span> <span class="mf">6.0</span><span class="o">*</span><span class="n">units</span><span class="o">.</span><span class="n">Ang</span>  <span class="c1"># Inner hysteretic radius for QM region</span>
<span class="n">qm_outer_radius</span> <span class="o">=</span> <span class="mf">8.0</span><span class="o">*</span><span class="n">units</span><span class="o">.</span><span class="n">Ang</span>  <span class="c1"># Inner hysteretic radius for QM region</span>
<span class="n">extrapolate_steps</span> <span class="o">=</span> <span class="mi">10</span>           <span class="c1"># Number of steps for predictor-corrector</span>
                                 <span class="c1"># interpolation and extrapolation</span>

<span class="c1"># ******* End of parameters *************</span>

<span class="n">set_fortran_indexing</span><span class="p">(</span><span class="kc">False</span><span class="p">)</span>

<span class="c1"># ********** Read input file ************</span>

<span class="nb">print</span> <span class="s1">&#39;Loading atoms from file </span><span class="si">%s</span><span class="s1">&#39;</span> <span class="o">%</span> <span class="n">input_file</span>
<span class="n">atoms</span> <span class="o">=</span> <span class="n">Atoms</span><span class="p">(</span><span class="n">input_file</span><span class="p">)</span>

<span class="n">orig_height</span> <span class="o">=</span> <span class="n">atoms</span><span class="o">.</span><span class="n">info</span><span class="p">[</span><span class="s1">&#39;OrigHeight&#39;</span><span class="p">]</span>
<span class="n">orig_crack_pos</span> <span class="o">=</span> <span class="n">atoms</span><span class="o">.</span><span class="n">info</span><span class="p">[</span><span class="s1">&#39;CrackPos&#39;</span><span class="p">]</span><span class="o">.</span><span class="n">copy</span><span class="p">()</span>

<span class="c1"># ***** Setup constraints *******</span>

<span class="n">top</span> <span class="o">=</span> <span class="n">atoms</span><span class="o">.</span><span class="n">positions</span><span class="p">[:,</span> <span class="mi">1</span><span class="p">]</span><span class="o">.</span><span class="n">max</span><span class="p">()</span>
<span class="n">bottom</span> <span class="o">=</span> <span class="n">atoms</span><span class="o">.</span><span class="n">positions</span><span class="p">[:,</span> <span class="mi">1</span><span class="p">]</span><span class="o">.</span><span class="n">min</span><span class="p">()</span>
<span class="n">left</span> <span class="o">=</span> <span class="n">atoms</span><span class="o">.</span><span class="n">positions</span><span class="p">[:,</span> <span class="mi">0</span><span class="p">]</span><span class="o">.</span><span class="n">min</span><span class="p">()</span>
<span class="n">right</span> <span class="o">=</span> <span class="n">atoms</span><span class="o">.</span><span class="n">positions</span><span class="p">[:,</span> <span class="mi">0</span><span class="p">]</span><span class="o">.</span><span class="n">max</span><span class="p">()</span>

<span class="c1"># fix atoms in the top and bottom rows</span>
<span class="n">fixed_mask</span> <span class="o">=</span> <span class="p">((</span><span class="nb">abs</span><span class="p">(</span><span class="n">atoms</span><span class="o">.</span><span class="n">positions</span><span class="p">[:,</span> <span class="mi">1</span><span class="p">]</span> <span class="o">-</span> <span class="n">top</span><span class="p">)</span> <span class="o">&lt;</span> <span class="mf">1.0</span><span class="p">)</span> <span class="o">|</span>
              <span class="p">(</span><span class="nb">abs</span><span class="p">(</span><span class="n">atoms</span><span class="o">.</span><span class="n">positions</span><span class="p">[:,</span> <span class="mi">1</span><span class="p">]</span> <span class="o">-</span> <span class="n">bottom</span><span class="p">)</span> <span class="o">&lt;</span> <span class="mf">1.0</span><span class="p">))</span>
<span class="n">fix_atoms</span> <span class="o">=</span> <span class="n">FixAtoms</span><span class="p">(</span><span class="n">mask</span><span class="o">=</span><span class="n">fixed_mask</span><span class="p">)</span>
<span class="nb">print</span><span class="p">(</span><span class="s1">&#39;Fixed </span><span class="si">%d</span><span class="s1"> atoms</span><span class="se">\n</span><span class="s1">&#39;</span> <span class="o">%</span> <span class="n">fixed_mask</span><span class="o">.</span><span class="n">sum</span><span class="p">())</span>
<span class="n">atoms</span><span class="o">.</span><span class="n">set_constraint</span><span class="p">([</span><span class="n">fix_atoms</span><span class="p">])</span>

<span class="c1"># Increase epsilon_yy applied to all atoms at constant strain rate</span>
<span class="n">strain_atoms</span> <span class="o">=</span> <span class="n">ConstantStrainRate</span><span class="p">(</span><span class="n">orig_height</span><span class="p">,</span> <span class="n">strain_rate</span><span class="o">*</span><span class="n">timestep</span><span class="p">)</span>

<span class="c1"># ******* Set up potentials and calculators ********</span>

<span class="n">mm_pot</span> <span class="o">=</span> <span class="n">Potential</span><span class="p">(</span><span class="n">mm_init_args</span><span class="p">,</span>
                   <span class="n">param_filename</span><span class="o">=</span><span class="n">param_file</span><span class="p">,</span>
                   <span class="n">cutoff_skin</span><span class="o">=</span><span class="n">cutoff_skin</span><span class="p">)</span>

<span class="c1"># Density functional tight binding (DFTB) potential</span>
<span class="n">qm_pot</span> <span class="o">=</span> <span class="n">Potential</span><span class="p">(</span><span class="n">qm_init_args</span><span class="p">,</span>
                   <span class="n">param_filename</span><span class="o">=</span><span class="n">param_file</span><span class="p">)</span>

<span class="c1"># Construct the QM/MM potential, which mixes QM and MM forces.</span>
<span class="c1"># The qm_args_str parameters control how the QM calculation is carried out:</span>
<span class="c1"># we use a single cluster, periodic in the z direction and terminated</span>
<span class="c1"># with hydrogen atoms. The positions of the outer layer of buffer atoms</span>
<span class="c1"># are not randomised.</span>
<span class="n">qmmm_pot</span> <span class="o">=</span> <span class="n">ForceMixingPotential</span><span class="p">(</span><span class="n">pot1</span><span class="o">=</span><span class="n">mm_pot</span><span class="p">,</span>
                                <span class="n">pot2</span><span class="o">=</span><span class="n">qm_pot</span><span class="p">,</span>
                                <span class="n">qm_args_str</span><span class="o">=</span><span class="s1">&#39;single_cluster cluster_periodic_z carve_cluster &#39;</span><span class="o">+</span>
                                            <span class="s1">&#39;terminate cluster_hopping=F randomise_buffer=F&#39;</span><span class="p">,</span>
                                <span class="n">fit_hops</span><span class="o">=</span><span class="mi">4</span><span class="p">,</span>
                                <span class="n">lotf_spring_hops</span><span class="o">=</span><span class="mi">3</span><span class="p">,</span>
                                <span class="n">hysteretic_buffer</span><span class="o">=</span><span class="kc">True</span><span class="p">,</span>
                                <span class="n">hysteretic_buffer_inner_radius</span><span class="o">=</span><span class="mf">7.0</span><span class="p">,</span>
                                <span class="n">hysteretic_buffer_outer_radius</span><span class="o">=</span><span class="mf">9.0</span><span class="p">,</span>
                                <span class="n">cluster_hopping_nneighb_only</span><span class="o">=</span><span class="kc">False</span><span class="p">,</span>
                                <span class="n">min_images_only</span><span class="o">=</span><span class="kc">True</span><span class="p">)</span>

<span class="c1"># Use the force mixing potential as the Atoms&#39; calculator</span>
<span class="n">atoms</span><span class="o">.</span><span class="n">set_calculator</span><span class="p">(</span><span class="n">qmmm_pot</span><span class="p">)</span>
<span class="n">qmmm_pot</span><span class="o">.</span><span class="n">atoms</span> <span class="o">=</span> <span class="n">atoms</span>

<span class="c1"># *** Set up the initial QM region ****</span>

<span class="n">qm_list</span> <span class="o">=</span> <span class="n">update_hysteretic_qm_region</span><span class="p">(</span><span class="n">atoms</span><span class="p">,</span> <span class="p">[],</span> <span class="n">orig_crack_pos</span><span class="p">,</span>
                                      <span class="n">qm_inner_radius</span><span class="p">,</span> <span class="n">qm_outer_radius</span><span class="p">)</span>

<span class="c1"># ********* Setup and run MD ***********</span>

<span class="c1"># Set the initial temperature to 2*simT: it will then equilibriate to</span>
<span class="c1"># simT, by the virial theorem</span>
<span class="n">MaxwellBoltzmannDistribution</span><span class="p">(</span><span class="n">atoms</span><span class="p">,</span> <span class="mf">2.0</span><span class="o">*</span><span class="n">sim_T</span><span class="p">)</span>

<span class="c1"># Initialise the dynamical system</span>
<span class="n">dynamics</span> <span class="o">=</span> <span class="n">LOTFDynamics</span><span class="p">(</span><span class="n">atoms</span><span class="p">,</span> <span class="n">timestep</span><span class="p">,</span> <span class="n">extrapolate_steps</span><span class="p">)</span>

<span class="c1"># array to store time averaged stress field</span>
<span class="n">avg_sigma</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">((</span><span class="nb">len</span><span class="p">(</span><span class="n">atoms</span><span class="p">),</span> <span class="mi">3</span><span class="p">,</span> <span class="mi">3</span><span class="p">))</span>

<span class="c1"># Print some information every time step</span>
<span class="k">def</span> <span class="nf">printstatus</span><span class="p">():</span>
    <span class="k">if</span> <span class="n">dynamics</span><span class="o">.</span><span class="n">nsteps</span> <span class="o">==</span> <span class="mi">1</span><span class="p">:</span>
        <span class="nb">print</span> <span class="s2">&quot;&quot;&quot;</span>
<span class="s2">State      Time/fs    Temp/K     Strain      G/(J/m^2)  CrackPos/A D(CrackPos)/A</span>
<span class="s2">---------------------------------------------------------------------------------&quot;&quot;&quot;</span>

    <span class="n">log_format</span> <span class="o">=</span> <span class="p">(</span><span class="s1">&#39;</span><span class="si">%(label)-4s%(time)12.1f%(temperature)12.6f</span><span class="s1">&#39;</span><span class="o">+</span>
                  <span class="s1">&#39;</span><span class="si">%(strain)12.5f%(G)12.4f%(crack_pos_x)12.2f</span><span class="s1">    (</span><span class="si">%(d_crack_pos_x)+5.2f</span><span class="s1">)&#39;</span><span class="p">)</span>

    <span class="n">atoms</span><span class="o">.</span><span class="n">info</span><span class="p">[</span><span class="s1">&#39;label&#39;</span><span class="p">]</span> <span class="o">=</span> <span class="n">dynamics</span><span class="o">.</span><span class="n">state_label</span>  <span class="c1"># Label for the status line</span>
    <span class="n">atoms</span><span class="o">.</span><span class="n">info</span><span class="p">[</span><span class="s1">&#39;time&#39;</span><span class="p">]</span> <span class="o">=</span> <span class="n">dynamics</span><span class="o">.</span><span class="n">get_time</span><span class="p">()</span><span class="o">/</span><span class="n">units</span><span class="o">.</span><span class="n">fs</span>
    <span class="n">atoms</span><span class="o">.</span><span class="n">info</span><span class="p">[</span><span class="s1">&#39;temperature&#39;</span><span class="p">]</span> <span class="o">=</span> <span class="p">(</span><span class="n">atoms</span><span class="o">.</span><span class="n">get_kinetic_energy</span><span class="p">()</span> <span class="o">/</span>
                                 <span class="p">(</span><span class="mf">1.5</span><span class="o">*</span><span class="n">units</span><span class="o">.</span><span class="n">kB</span><span class="o">*</span><span class="nb">len</span><span class="p">(</span><span class="n">atoms</span><span class="p">)))</span>
    <span class="n">atoms</span><span class="o">.</span><span class="n">info</span><span class="p">[</span><span class="s1">&#39;strain&#39;</span><span class="p">]</span> <span class="o">=</span> <span class="n">get_strain</span><span class="p">(</span><span class="n">atoms</span><span class="p">)</span>
    <span class="n">atoms</span><span class="o">.</span><span class="n">info</span><span class="p">[</span><span class="s1">&#39;G&#39;</span><span class="p">]</span> <span class="o">=</span> <span class="n">get_energy_release_rate</span><span class="p">(</span><span class="n">atoms</span><span class="p">)</span><span class="o">/</span><span class="p">(</span><span class="n">units</span><span class="o">.</span><span class="n">J</span><span class="o">/</span><span class="n">units</span><span class="o">.</span><span class="n">m</span><span class="o">**</span><span class="mi">2</span><span class="p">)</span>

    <span class="n">crack_pos</span> <span class="o">=</span> <span class="n">find_crack_tip_stress_field</span><span class="p">(</span><span class="n">atoms</span><span class="p">,</span> <span class="n">calc</span><span class="o">=</span><span class="n">mm_pot</span><span class="p">,</span>
                                            <span class="n">avg_sigma</span><span class="o">=</span><span class="n">avg_sigma</span><span class="p">)</span>
    <span class="n">atoms</span><span class="o">.</span><span class="n">info</span><span class="p">[</span><span class="s1">&#39;crack_pos_x&#39;</span><span class="p">]</span> <span class="o">=</span> <span class="n">crack_pos</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span>
    <span class="n">atoms</span><span class="o">.</span><span class="n">info</span><span class="p">[</span><span class="s1">&#39;d_crack_pos_x&#39;</span><span class="p">]</span> <span class="o">=</span> <span class="n">crack_pos</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span> <span class="o">-</span> <span class="n">orig_crack_pos</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span>

    <span class="nb">print</span> <span class="n">log_format</span> <span class="o">%</span> <span class="n">atoms</span><span class="o">.</span><span class="n">info</span>

<span class="n">dynamics</span><span class="o">.</span><span class="n">attach</span><span class="p">(</span><span class="n">printstatus</span><span class="p">)</span>

<span class="k">def</span> <span class="nf">atom_straining</span><span class="p">(</span><span class="n">atoms</span><span class="p">):</span>
    <span class="n">crack_pos</span> <span class="o">=</span> <span class="n">find_crack_tip_stress_field</span><span class="p">(</span><span class="n">atoms</span><span class="p">,</span> <span class="n">calc</span><span class="o">=</span><span class="n">mm_pot</span><span class="p">,</span>
                                            <span class="n">avg_sigma</span><span class="o">=</span><span class="n">avg_sigma</span><span class="p">)</span>
    <span class="c1"># keep straining until the crack tip has advanced to tip_move_tol</span>
    <span class="k">if</span> <span class="ow">not</span> <span class="n">atoms</span><span class="o">.</span><span class="n">info</span><span class="p">[</span><span class="s1">&#39;is_cracked&#39;</span><span class="p">]</span> <span class="ow">and</span> <span class="p">(</span><span class="n">crack_pos</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span> <span class="o">-</span> <span class="n">orig_crack_pos</span><span class="p">[</span><span class="mi">0</span><span class="p">])</span> <span class="o">&lt;</span> <span class="n">tip_move_tol</span><span class="p">:</span>
      <span class="n">strain_atoms</span><span class="o">.</span><span class="n">apply_strain</span><span class="p">(</span><span class="n">atoms</span><span class="p">)</span>
    <span class="k">elif</span> <span class="ow">not</span> <span class="n">atoms</span><span class="o">.</span><span class="n">info</span><span class="p">[</span><span class="s1">&#39;is_cracked&#39;</span><span class="p">]:</span>
      <span class="n">atoms</span><span class="o">.</span><span class="n">info</span><span class="p">[</span><span class="s1">&#39;is_cracked&#39;</span><span class="p">]</span> <span class="o">=</span> <span class="kc">True</span>

<span class="n">dynamics</span><span class="o">.</span><span class="n">attach</span><span class="p">(</span><span class="n">atom_straining</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="n">dynamics</span><span class="o">.</span><span class="n">atoms</span><span class="p">)</span>

<span class="c1"># Function to update the QM region at the beginning of each extrapolation cycle</span>
<span class="k">def</span> <span class="nf">update_qm_region</span><span class="p">(</span><span class="n">atoms</span><span class="p">):</span>
   <span class="n">crack_pos</span> <span class="o">=</span> <span class="n">find_crack_tip_stress_field</span><span class="p">(</span><span class="n">atoms</span><span class="p">,</span> <span class="n">calc</span><span class="o">=</span><span class="n">mm_pot</span><span class="p">,</span>
                                           <span class="n">avg_sigma</span><span class="o">=</span><span class="n">avg_sigma</span><span class="p">)</span>
   <span class="n">qm_list</span> <span class="o">=</span> <span class="n">qmmm_pot</span><span class="o">.</span><span class="n">get_qm_atoms</span><span class="p">(</span><span class="n">atoms</span><span class="p">)</span>
   <span class="n">qm_list</span> <span class="o">=</span> <span class="n">update_hysteretic_qm_region</span><span class="p">(</span><span class="n">atoms</span><span class="p">,</span> <span class="n">qm_list</span><span class="p">,</span> <span class="n">crack_pos</span><span class="p">,</span>
                                         <span class="n">qm_inner_radius</span><span class="p">,</span> <span class="n">qm_outer_radius</span><span class="p">)</span>
   <span class="n">qmmm_pot</span><span class="o">.</span><span class="n">set_qm_atoms</span><span class="p">(</span><span class="n">qm_list</span><span class="p">,</span> <span class="n">atoms</span><span class="p">)</span>
   <span class="c1">#assert (atoms.hybrid == 1).sum() == len(qm_list)</span>

<span class="n">dynamics</span><span class="o">.</span><span class="n">set_qm_update_func</span><span class="p">(</span><span class="n">update_qm_region</span><span class="p">)</span>


<span class="c1"># Save frames to the trajectory every `traj_interval` time steps</span>
<span class="c1"># but only when interpolating</span>
<span class="n">trajectory</span> <span class="o">=</span> <span class="n">AtomsWriter</span><span class="p">(</span><span class="n">traj_file</span><span class="p">)</span>

<span class="k">def</span> <span class="nf">traj_writer</span><span class="p">(</span><span class="n">dynamics</span><span class="p">):</span>
   <span class="k">if</span> <span class="n">dynamics</span><span class="o">.</span><span class="n">state</span> <span class="o">==</span> <span class="n">LOTFDynamics</span><span class="o">.</span><span class="n">Interpolation</span><span class="p">:</span>
      <span class="c1"># copy time-averaged stress into Atoms so it will be written to trajectory file</span>
      <span class="n">dynamics</span><span class="o">.</span><span class="n">atoms</span><span class="o">.</span><span class="n">set_array</span><span class="p">(</span><span class="s1">&#39;avg_sigma&#39;</span><span class="p">,</span> <span class="n">avg_sigma</span><span class="o">.</span><span class="n">reshape</span><span class="p">((</span><span class="nb">len</span><span class="p">(</span><span class="n">atoms</span><span class="p">),</span> <span class="mi">9</span><span class="p">)))</span>
      <span class="n">trajectory</span><span class="o">.</span><span class="n">write</span><span class="p">(</span><span class="n">dynamics</span><span class="o">.</span><span class="n">atoms</span><span class="p">)</span>

<span class="n">dynamics</span><span class="o">.</span><span class="n">attach</span><span class="p">(</span><span class="n">traj_writer</span><span class="p">,</span> <span class="n">traj_interval</span><span class="p">,</span> <span class="n">dynamics</span><span class="p">)</span>

<span class="c1"># Start running!</span>
<span class="n">dynamics</span><span class="o">.</span><span class="n">run</span><span class="p">(</span><span class="n">nsteps</span><span class="p">)</span>
</pre></div>
</div>
</div>
</div>


           </div>
           
          </div>
          <footer>
  
    <div class="rst-footer-buttons" role="navigation" aria-label="footer navigation">
      
        <a href="adaptive-qmmm-references.html" class="btn btn-neutral float-right" title="References" accesskey="n" rel="next">Next <span class="fa fa-arrow-circle-right"></span></a>
      
      
        <a href="adaptive-qmmm-step3.html" class="btn btn-neutral float-left" title="Step 3: LOTF hybrid MD simulation of fracture in Si" accesskey="p" rel="prev"><span class="fa fa-arrow-circle-left"></span> Previous</a>
      
    </div>
  

  <hr/>

  <div role="contentinfo">
    <p>
        &copy; Copyright 2008-2016, James Kermode

    </p>
  </div>
  Built with <a href="http://sphinx-doc.org/">Sphinx</a> using a <a href="https://github.com/rtfd/sphinx_rtd_theme">theme</a> provided by <a href="https://readthedocs.org">Read the Docs</a>. 

</footer>

        </div>
      </div>

    </section>

  </div>
  


  <script type="text/javascript">
      jQuery(function () {
          SphinxRtdTheme.Navigation.enable(true);
      });
  </script>

  
  
    
   

</body>
</html>